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We investigate the effect that the average backreaction of structure formation has on the dynamics 
of the cosmological expansion, within the concordance model. Our approach in the Poisson gauge is 
fully consistent up to second-order in a perturbative expansion about a flat Friedmann background, 
including a cosmological constant. We discuss the key length scales which are inherent in any 
averaging procedure of this kind. We identify an intrinsic homogeneity scale that arises from the 
averaging procedure, beyond which a residual offset remains in the expansion rate and deceleration 
parameter. In the case of the deceleration parameter, this can lead to a quite large increase in the 
value, and may therefore have important ramifications for dark energy measurements, even if the 
underlying nature of dark energy is a cosmological constant. We give the intrinsic variance that 
affects the value of the effective Hubble rate and deceleration parameter. These considerations serve 
to add extra intrinsic errors to our determination of the cosmological parameters, and, in particular, 
may render attempts to measure the Hubble constant to percent precision overly optimistic. 
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I. INTRODUCTION 

The universe appears to be close to homogeneous and isotropic, on average, on large scales, but it exhibits a very 
clumpy distribution of matter on small scales. To account for this structure, the standard cosmological model relics 
on the separation of the geometry of spacetime into a perfectly homogeneous and isotropic Friedmann-Lemaitre- 
Robcrtson- Walker (FLRW) background describing the large scale properties of the universe, such as the expansion 
rate, and small fluctuations around this background solution. This provides a straightforward perturbative treatment 
of the growth of structure under the influence of gravitation. The explicit construction of the background by a 
smoothing or averaging procedure applied to the clumpy Universe is often ignored, and the background appears as 
an artificial mathematical object used to perform the calculations of gauge invariant quantities characterising the 
physical properties of the clumpy universe. 

The essence of this 'averaging problem' comes when we try to match the late time universe today, which is full of 
structure, to the early time universe, which isn't. At the end of inflation we are left with a universe with curvature 
characterised by some constant fci n f (= 0,±1 in some units), and cosmological constant, Aj n f, which are fixed for all 
time (and might be zero), and perturbations which are of tiny amplitude and well outside the Hubble radius; there is 
no averaging problem at this time, and the idea of background plus perturbations is very natural and simple to define. 
Fast-forward to today, where structures are non-linear, are inside the Hubble radius, and many have broken away from 
the cosmic expansion altogether. We may still apparently describe the universe as FLRW plus perturbations to high 
accuracy; that is, it is natural and seemingly correct to define a FLRW background, but it is implicitly assumed that 
this background is the same one that we are left with at the end of inflation, in terms of k- m { and A; n f. Mathematically 
we can follow a model from inflation to today, but when we try to fit our models to observations to describe our 
local universe we arc implicitly smoothing over structure, and this can contaminate what we think our inflationary 
background FLRW model should be. Indeed, it is not clear that the background smoothed model should actually 
obey the field equations at all. Within the standard paradigm, then, the averaging problem also becomes a fitting 
problem; are the background parameters we are fitting with the CMB actually the same as those when fitting SNIa? 
(See [1] for a discussion of these issues.) 

Because of the non-linearity of the Einstein field equations, the explicit construction of a homogeneous background 
is far from trivial and it has been know for a long time that the local fluctuations may influence the way the Universe 
behaves on average [2]; this effect is usually dubbed backreaction and has started to be investigated in detail (see, 
e.g., [3-26] and references therein), often as a possible solution to the dark energy problem itself. The problem 
is mathematical: if we have an inhomogeneous matter distribution in some spacetime, and we try to calculate a 
homogeneous 'background' by smoothing the matter content and calculating the new smoothed metric, we get a 
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different answer than if we smooth the metric directly; this difference is usually termed backreaction, in this context. 

In this work, we explore the changes to the background due to the presence of perturbations by explicitly calculating 
the effects of backreaction up to second-order in perturbation theory. Similar investigations have been done in the 
past in the synchronous gauge [17-19], and in the Poisson gauge (see e.g. [5, 6, 20, 21, 23]). The authors of these 
previous works have mainly considered only terms which are quadratic in linear quantities (ignoring the second-order 
Bardeen potentials), and/or a pure Einstein de-Sitter Universe (ignoring the cosmological constant). 

Here, we present an analysis which is consistent up to second order, both for Einstein de-Sitter and the concordance 
model - flat ACDM - as underlying background solutions. To perform our averages we argue that the rest frame of 
the gravitational field is a natural choice of spatial hypersurface, and that this may be unambiguously defined by the 
vanishing of the magnetic part of the Weyl tensor, where this is possible. We then define the effective expansion rate 
through the average of the matter fluid expansion projected onto these spatial hypersurfaces, rather than through 
the expansion 4-scalar for matter (as in [5, 6]) or through the expansion of the 'observer' (defined below) flow lines 
(as in [20-23] - though as we show below, this is in fact the rest frame of the gravitational potential). Our choice 
can be justified by the fact that the averaging process used in this work is frame (observer) dependent, so that, for 
a matter of consistency, one wants to retain the observer dependence on the quantities that are averaged also. On 
the one hand, the usual expansion 4-scalar for matter is the observed expansion only for observers comoving with the 
fluid so that it is only appropriate for averaging with respect to a set of comoving observers. On the other hand, the 
expansion of the observer flow lines is clearly frame dependent, but it does not encompass any information about the 
matter flow, apart from its background part; the averaged homogeneity being a characteristic of the total matter flow, 
it seems then more natural to retain a quantity that takes into account the fluctuations in the matter distribution 
and its peculiar velocity with respect to the observer's rest-frame. 

Moreover, we include an important characteristic effect of backreaction, viz. the existence of an intrinsic variance 
on quantities such that the Hubble rate and the deceleration parameter (as noted by [6, 18]). This is important 
because it gives us our '1 — a error' which is an intrinsic unknown when identifying a background based on Gaussian 
perturbations around FLRW. 

The paper is organized as follows. In Section 2 we describe the averaging formalism used in the rest of the paper and 
based on a recent generalization [26] of the standard averaging procedure (defined in a comoving coordinate system) 
to arbitrary coordinate systems. In Section 3, we apply this formalism to the averaging of the cosmological model 
in the Poisson gauge up to second order. In Section 4, we present the effects on the Hubble rate and deceleration 
parameter, with a particular emphasis on the importance of the different scales involved (smoothing and averaging) 
in the process. We find that the effect on the quantities themselves is small (as expected in perturbation theory) 
but that it can be quite large in terms of the variance affecting these quantities. Our results are quantitatively in 
agreement with the results obtained in the synchronous gauge [17-19]. Finally, Section 5 summarizes the results and 
addresses quickly possible future devclopements of this work. 

II. AVERAGING FORMALISM IN AN ARBITRARY COORDINATE SYSTEM 

In this paper, we are concerned with estimating the backreaction effect in the Poisson gauge of perturbation theory. 
In this gauge, the 4- velocity of the matter fluid is tilted with respect to the timelike normal vector of the coordinate 
system, and this makes it difficult to use the standard averaging procedure (see [3]), that has been developed with 
respect to observers comoving with the cosmic matter fluid. A recent work [26] has generalized the averaging procedure 
to an arbitrary coordinate system, which allows us to estimate consistently the backreaction effect in the Poisson gauge. 
This section briefly presents this formalism and discusses the assumptions at the root of it. Essentially, we introduce 
two velocity fields, one in the rest-frame of the matter content, and another arbitrary one on which we fix an arbitrary 
family of 'observers'. The role of the observers is to introduce spacelike hypersurfaces on which we perform our 
average; these hypersurfaces are tilted with respect to the spacelike hypersurfaces defined by the fluid (tilted in the 
sense that their normal vectors are tilted). We fix the coordinates with respect to the observers, and so loosely refer 
to this as the coordinate frame. 

Throughout the paper, we will suppose that gravitation is well described by general relativity on all scales and that 
the cosmic matter fluid can be considered as a perfect fluid. Moreover, Latin letters of the beginning of the alphabet 
(a, b,c,...,h) will denote spacetime indices, and Latin letters in the middle the alphabet (i, j,k, ...) will denote spatial 
indices. 
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A. 1 + 3 threading of spacetime 

We consider a set of observers 0(jp) denned at each point of the spacetime manifold p £ A4, and characterized 
by a unit 4- velocity field n a that is everywhere timclikc and future directed, i.e. n a n a = — 1, with zero vorticity. 
This 4-velocity field induces a natural foliation of spacetime by a continuous set of space-like hypersurfaces locally 
orthogonal to n a . Then we can define the projection tensor field onto these hypersurfaces as h a b = g a b + n a rib that is 
a well-defined Ricmannian metric for these hypersurfaces. The line clement can then be written, with respect to this 
foliation: 

ds 2 = -{N 2 - NiN l )dt 2 + 2N,dtdx l + h^dx'dx 1 , (1) 

where we have introduced respectively the lapse function N(x a ) and the shift 3- vector N l (x a ) such that the components 
of the 4-velocity read: 

n a = , n a = N(-l, 0,0,0) . (2) 

For the purposes of this paper, the hypersurfaces orthogonal to n a are characterized by two quantities: 

• the intrinsic curvature of the hypersurfaces: 1Z = h ab lZ a b, where 1Z a b is the 3-Ricci curvature of the hypersurfaces; 

• the extrinsic curvature (or second fundamental form): K a b = —h a h^n c -d that encodes the way the hypersurfaces 
are embedded in the manifold M. 

In the following, we will assume that the matter content can be well described by a perfect fluid (not necessarily 
irrotational) of energy density p(x a ), pressure p(x a ) and 4-velocity u a (x b ) (with u a u a — —1), so that its stress-energy 
tensor reads: 

Tab ={p + P)u a u b + PQab ■ (3) 

Note that in this work, the 4-velocity of matter u a is not necessarily aligned with the 4-velocity of the observers n a , 
so that there exists a vector field v a corresponding to the relative velocity of the matter fluid with respect to the 
fundamental observers. v a is space-like and orthogonal to n a (v a n a — 0) and one has: 

u a = 7 (n a + v a ) with 7 = 1 , (4) 

where 7 is the usual Lorentz factor and v 2 — v a v a . Thanks to the 1+3 foliation, the Einstein field equations can 
be separated in two different sets: the constraint equations that have to be satisfied on every hypersurface, and the 
evolution equations, that prescribe how the fields (h a b,K a b) evolve from one hypersurface to another infinitesimally 
close. The Hamiltonian constraint reads: 

TZ - K)K{ +K 2 = 167rGe + 2A , where e = T ab n a n b = j 2 p + ( 7 2 - l)p , (5) 

where K = K\. The momentum constraint is: 

WiK) - VjK = 8wGJj , where J, = -T ab n a h) = 7 2 (p + p) Vj , (6) 

where we have defined the projected covariant 3-derivative on the spatial hypersurfaces of any tensor field t a - c d f -: 
^ 'dt°"" c a ...c = ^d^a' ■■■^c'^ l a ' ■■■^■s ' ^ et a ' ' c 'a>' c " ■ The evolution equation for the first fundamental form reads: 

jjdthij = -2Kij + ^VyJVfl , (7) 
and for the second fundamental form, one has: 

^8 t K) = TZ) + KK) - AS} - ^VjV'N + 1 (K^ 3 N k - K*\7 k N* + N k \7 k K<) - 8nG ^ + * (e - S k k )6^j , (8) 

with Sij = j 2 pviVj+p(hij + j 2 ViVj). These equations have to be supplemented by the energy-momentum conservation 
for the matter fluid: V a T fc a = 0. We will now introduce the standard decomposition for the covariant spatial derivatives 
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of the 4-vectors in terms of their trace, symmetric trace-free and antisymmetric parts. Writing / = n a V a f for any 
quantity /, one has: 

\7 a n b = -n a h b + ^hab + Sab (9) 

with £ ee h c a h da \7 c n d and £ ao ee h c a hfV {c n d) - ^h ab ; 
V a u h = -n a (2n b + vvn b ) - -y\7 a vn b + ^9h ab + a ab + uj ab (10) 

with 6 ee h c a h da \7 c u d , <j ab = h a h d b ^ {c u d) - ^9h ab and iv ab = h a h d b ^y c u d] ; 

\7 a v b = -n a (v b + X b ) + Y a n b +^Kh ab + f3 ab + W ab (11) 

with k ee h c a h da \7 c v d , [3 ab ee h c a hf\7 (c v d) - ^nh ab 

and X a ee n b h c a V c v b , Y b = n c h a b XI c v a , W ab = h c a hfV [c v d] . 

In these relations, £, and k denote the isotropic expansion rates of the 4- velocities n a , u a and of the peculiar velocity 
v a respectively, and T, ab , a ab and (3 ab their shears, with respect to the threading of spacetime induced by the vector 
n a . oj ab and W ab are, respectively, the vorticities 1 of u a and v a in this same foliation. These quantities are those 
measured by the observers with 4- velocities n a in their instantaneous rest-frame. In particular 9, a ab and uo ab differ 
from the usual expansion, shear and vorticity of the matter fluid as measured by observers comoving with this matter 
fluid (by acceleration terms essentially), that are defined by the decomposition of (g ac + u a u c ){g bd + u h u d )V c u d . For 
example, the expansions are linked by the relation: 

9 ee V a u a = 9 + j{j 2 v a v a - n a v a ) . (12) 

Using (4), one can relate these quantities as follows: 

£ = l~ 1 9-n- 1 2 B (13) 

S Q 6 = 7 _1(T ab - /?ab - 7 2 (yB ( ab ) - -Bh ab 

W ab = 7~ 1 w ab - -f 2 B [ab] , (15) 
where we have introduced the tensor: 

B ab = ^K(v a n b + v a v b ) + f3 ca v c n b + (3 ca v c v b + W ca v c n b + W ca v c v b , (16) 

whose trace is given by B — \nv 2 + /3 ab v a v b . In our notation, angular and round brackets denote the antisymmetric 
and symmetric parts, respectively, of a tensor projected with h ab - Let's finally introduce the following notation for 
convenience: 

B = -7K-7 3 £ (17) 



(14) 



UBij = -iPa - 7 3 (B{ij) - \ Bh i 



(18) 



so that: 



Z = 1 - 1 {9 + 9 B ), (19) 

=1~ X {<rn+<TBij) ■ (20) 

The non-local, free gravitational field is described by the Weyl tensor. Given a timelike vector this is split into 
electric and magnetic parts. For example, with respect to n a these are 

= C acbd n c n d and H<£> = *C acbd n c n d , (21) 



1 Since it defines the foliation, n a is vorticity free by definition. 
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where C a bcd is the Weyl tensor and *C a bcd is its dual. Analogous definitions exist for the vector field u a . This means 
that observers in the frame of the fluid and observers in the coordinate frame observe this electric-magnetc split 
differently (see [27] for the transformation relations between the two), analogously to boosted observers measuring 
different electric and magnetic parts of the electromagnetic field. In particular, in certain gravitational fields there 
may exist a special frame whereby one of these two components vanishes. For example, in so-called silent universes 
which arc not conformally flat, there exists a preferred frame in which the magnetic part of the Weyl tensor is zero - 
such a frame may be considered the rest-frame of the gravitational field. In spacetimes where this is possible, it is 
unique as follows from the transformation laws in [27], and there exist (at least) two physical, well motivated, frames: 
the rest-frame of the fluid, and the rest-frame of the non-local gravitational field. We return to this below. 

B. General averaging procedure 

The non-locality of the spatial average that is usually calculated manifests itself problematically when interpreting 
averaged quantities. What do they mean? In which spacetime do they exist - the rough or the smooth? - representing 
physical things (i.e., where are they tcnsorial objects)? It is tempting to average the fluid expansion of a lumpy 
spacetime, say, and to think of this as actually being in some sense the 'averaged fluid expansion rate'. Normally, of 
course, the expansion rate of a fluid is a covariantly defined local object, and so unambiguous when defined in the 
local rest-frame of the fluid. But it is not understood how to define the rest-frame of the non-local smoothed fluid in 
a covariant way, and the 'average expansion rate' picks up this ambiguity. 

Furthermore, any definition of an averaged expansion rate in the spacetime in which we started is not covariant 
from a 4-dimensional perspective because the average is with respect to spatial hypersurfaces defined by n a in the 
un-smoothed spacetime, and so implicitly rely on some coordinates, and, hence, a mapping between the unsmoothcd 
and smoothed spacetimes. How do we choose these coordinates? Scalar averaging approaches give us access to some 
averaged quantities and backreaction terms, but not the spacetime in which they exist as the objects they are supposed 
to represent. Analogously to gauge freedom, this ambiguity leaves implicit choices to make about what objects we 
consider important, as well as what they mean. Previous analyses have fixed these freedoms in one way or another. 

Ideally, we would like to construct a smooth FLRW 'spacetime' from an lumpy inhomogeneous spacetime by 
averaging over structure. This would, in principle, have a metric 

ds 2 ff = -dr 2 + o^-dy'dy* , (22) 

where r is the cosmic time and a-p (t) a scale factor, the subscript T> indicating that it has been obtained at a certain 
spatial scale characteristic of a compact spatial domain V, which is large enough so that a homogeneity scale has been 
reached; in this case 7^ will be a metric of constant curvature. This is not to imply that this would be a spacetime in 
the usual sense because the normal observational relations may not follow directly from this FLRW metric: these may 
have to be calculated separately [28, 29]. However, in the context of perturbation theory the effect of renormalisation 
of the background appears naturally as a first step in this procedure. 

In the context of averaging a perturbed spacetime we consider below, we can imagine choosing coordinates which 
straddle both the rough and smooth spacetimes. In particular, we can choose our coordinates in the rough spacetime 
such that they become the ones we want in the smooth spacetime. We will choose our coordinates such that the time 
coordinate in the rough spacetime becomes the proper time coordinate in the smoothed one: that is, we will set r = t 
as well as y % — x*. (We will also set iVj = 0.) This effectively de-synchronises the clocks between the rough (with 
proper time J Ndt when Ni = 0) and smoothed (proper time t) spacetimes. The averaging operator we define is 
simply the Riemannian average over the domain V in the surfaces orthogonal to n a (i.e., t =const.): 

(iP) v = ±-j Jd 3 x^(t,x*) , (23) 

and is well defined for any scalar function ip. This choice has the property that the commutator between partial time 
derivatives and spatial averages which reads, when N l = 0, 

[dr, x*) = (N^) v - (N0 V W v , (24) 

is zero when we consider perturbed Einstcin-dc Sitter models below at second-order (that is, the expansion of the 
spatial hypersurfaces £ when scaled by N is, for the perturbed model we consider below, | N£ = H — <f> — 2<f><I> — j^ 2 ) ; 
with 4> = as it is for pure dust, the commutator vanishes at second-order). 

Of the three expansion rates we have introduced, £ tells us the expansion of the coordinate grid, and so is not 
physically attached to the fluid. measures the fluid expansion rate in its own rest frame, and is the sensible choice 
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of expansion rate with which to characterise the rough spacetimc. However, as we have mentioned, after smoothing, 
the rest-frame of the fluid will change in a way which is not yet known, and will depend on the domain. Any 
expansion rate we try to investigate must take this into account and so allow for a tilt between the fluid and the 
normal to the hypersurfaces we use to average. The expansion rate 9, which is the trace of the expansion tensor of 
the fluid projected into the coordinate rest-frame in which the averaging takes place is the most natural choice when 
accounting for peculiar velocities in this way. We shall define our observers, which define the spatial hypersurfaces on 
which the averaging takes place, by the rest frame of the gravitational field; that is, the frame in which — 0. At 
second-order in a perturbation expansion, this is different from the rest-frame of the matter. 
Thus, in the rough spacetime, when we consider the length-scale I associated with 9, we have 

!„ n ^ , „ 1 d£ 1 dl 

Hence, if t represents the proper time in the smoothed spacetime, we can define a pre-synchronised, smoothed, Hubble 
parameter using N9 as [4] 



H v = l -{N6) V = -L J Jd 3 xN6 



(25) 



, d?ay 
a-D 



We may think of this as the average Hubble parameter which preserves the length-scale £ after smoothing, according 
to the pre-chosen proper time in the smoothed spacetime. This is referred to as the scaled t-Hubble parameter in [4], 
where it is introduced to preserve the structure of the averaged equations. This is not a unique choice and we refer 
to [10] for a detailed discussion. In particular, this Hubble parameter changes under a re-scaling of t; in the context 
of perturbations of FLRW, we use t to be the proper time in the background - if we were to use conformal time, say, 
Ht> would be an averaged conformal Hubble parameter. We can use this to then define the effective scale factor for 
the averaged model as the function a-p(t) obeying: 

H V (26) 
a v 

Using the commutation relation, one can average the scalar part of the Einstein field equations to obtain a set of two 
equations giving the behavior of the effective scale factor ax>{t) (see [26] with TVj = 0): 

ml = 16ttG ((l i N 2 p) v + ( 7 2 ( 7 2 - l)N 2 p) v ) + 2A (N 2 j 2 ) v - (j 2 N 2 TZ) v - Q v + C v (27) 

-4nG (N 2 j ( 7 2 p(l + v 2 ) + 3p(l + 2 7 V))) I , + A (N 2 1 ) v 
+ Qv + Vv+^ V + T v - C v , (28) 
where we have defined the standard kinematical backreaction: 

Q v ee I (((N9) 2 ) v (N9) 2 V ) 2 (N 2 a 2 ) v , (29) 
and additional backreaction terms as: 

C V ee 2 (N 2 a 2 B ) v I ((N9 B ) 2 ) V \ (N 2 99 B ) V (30) 

V v = (9d t N) v + ^NV k V k N^ (31) 

IC V = (N 2 1 - 1 9 B 9) V -3(N 1 - 1 9 B ) V H V 

- (N8 t 9 B ) v + (N 2 ^e) v + {N 2 1 - 1 j9 B ) v - 2 (N 2 9 2 B ) V (32) 

Fv - l(N 2 9 2 ( 1 - 1 -l)) v -2(N 2 a 2 ( 1 - 1 -l)) v -(N9) v (N9( 1 - 1 -l)) v 

-\ (N 2 9 2 B (^ - l)) v - \ (N 2 66 B ( 7 -i l)) v 2 (N 2 a%(^ l)) v . (33) 



III. AVERAGING PERTURBED FLRW MODELS 



We consider the backreaction effect in a perturbed FLRW model, with a flat background and a cosmological constant. 
We are interested in the backreaction effect at late times, and so assume the matter in our model is comoving cold dark 
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matter plus baryons. We shall consider scalar modes up to second-order, and ignore vectors and tensors throughout, 
as the first-order contributions are small and it has been shown that the vectors and tensors induced by first-order 
scalars are also sub-dominant, although they could provide a slight correction to the results presented here [30-34]. 
In the Poisson gauge 2 the metric reads [35] 

ds 2 = - (l + 2$ + * (2) ) dt 2 + a 2 (l - 24- - *( 2 >) (SydscW. (34) 

The first-order scalar perturbations are given by < i ) , v l', and the second-order by & 2 \ ty( 2 \ In this form we have the 
metric in its Newtonian form, which we may think of as the local rest-frame of the gravitational field. Comparing 
with the general metric given above, we are in a gauge with zero shift vector, and TV 2 = (l + 2<f> + <f>( 2 )) . Moreover, 
the matter fluid has a peculiar velocity v a — (0,v l )/a where vi — ^di(2v^ + v^), where is the first order part 
and the second order part of the velocity potential. The spatial metric is obviously — a 2 (l — 2^ — "ft 2 )) 5ij. 
The Poisson gauge is particularly elegant for scalar perturbations because with n a defined orthogonal to the spatial 
metric hij, the Weyl tensor becomes 



E[f = \ {h i a h ] b - ^hijh^ jv„V& $ + *-$ 2 -* 2 + i(V 2 ) + *( 2 )) +V„$V6$-V„*V6*| 



(35) 



- 0. (36) 

In the rest frame n a , then, the gravitational field is silent, and, with ^ = <f> is a pure potential field. Hence, n a may 
be considered as the rest-frame of the gravitational field, and so defines natural hypersurfaces with which to perform 
our averages. By contrast, in the frame u a the Weyl tensor has non-zero H a b [27]. 

As time coordinates we shall use conformal time, r\, proper time in the background, t, and background redshift 
z = l/a — 1, related by adrj = dt = —dz/H(l + z), interchangeably, where the background Hubble rate H = a/a = 
l/a(da/dt) is given by 

H(z) 2 = H 2 [fi (l + z) 3 + 1 - fi ] (37) 

and f2 = £l m (z = 0) is the present-day matter density parameter. Note that in the perturbed spacetime the parameter 
z is just a time coordinate, even though we refer to it as redshift. The Raychaudhuri equation in the background, 
H = — |i? 2 fi m , gives the deceleration parameter 

. . la „ 1 + z dH 3_ , . 

q ~ i{z) = -JPa=- l+ H{z)^ = - l+ 2 nm{2) > (38) 

where 

to m (z) = "° (1 + Z)3 ^, (39) 

[«o(l + Z ) 3 + l-^o] 1/2 

gives the evolution of the density parameter; the density parameter associated with the cosmological constant is 
^a(z) = 1 - fi m (z). 

For a single fluid with zero pressure and no anisotropic stress $ = <&, and <f> obeys the Bardeen equation 

+ m& + a 2 A$ = = $ + 4H$ + A$ . (40) 

and ' = d/dij, and H = a' /a is the conformal Hubble rate. All first-order quantities can be derived from <J>; for 
example, 

w(1) --3^(* + ^)' ^ 

and $ is the source of the second-order scalars. The solution to the growing mode of the Bardeen equation may be 
written as 

$(»7,x)=0(T7)$o(aO (42) 



2 In this work, we call a gauge the choice of both the frame, i.e. n a , and the coordinate set on hypersurfaces orthogonal to n a . 
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where &o( x ) is the Bardeen potential today (rj = rjo, z = 0) and g(rj) is the growth suppression factor, which may be 
approximated, in terms of rcdshift, as [36, 37] 



1 + 2^™( Z ) 



(43) 



and goo is chosen so that g(z = 0) = 1. 

We define our Fourier transform as (suppressing any temporal quantities) 



*(aj) 



(2tt) 3 / 2 

where $*(fe) = <&(— fe). The power spectrum of <& is defined by 

2ir 2 



d 3 k<I>(k)e 



ikx 



$(fc)$(fc') = ^Vt>{k)5{k + k'), 



(44) 



(45) 



where an overbar denotes an ensemble average. Assuming scale-invariant initial conditions from inflation, this is given 

by 



V^{z,k) 



3A R 

5goo 



5 (*) 2 T(fc) 2 



(46) 



where T{k) is the normalised transfer function, A 2 ^ is the primordial power of the curvature perturbation, with [50] 
A| w 2.41 x 1CT 9 at a scale k CMB = 0.002Mpc _1 . 

The second-order solutions for and <I>( 2 ) are given by [35]. We quote their results directly: 



* (2) (?7,a0 



+■ ( ^2(1?) - -g(v)9r, 
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(a n i - i)g{v)g m H 



v- 2 (d l $o^o) - 3v- 4 a 4 ^ (9^0^*0) 



$ (2) (t?,:e) 



Si(r?) + 4 5 2 (?7) - 2g(r ] )g m - — (a n] - 1)3(77)3™ ) <I>n 



(47) 



•82(7?) + o5 2 (?7) e(7?) + - - -3(77)3 



V- 2 (&$ Q d& ) - (d i $ d J $o) 



+B 3 ( v )V- 2 d l d^d^ Q d^o) + B 4 (77)9^oft$o : 



(48) 



where Bi(rj) —H 2 (Jo + 3fio/2) 1 Bi(i]) with the following definitions 



Bi(v) 



df}H 2 (v)(m-l) 2 C(v,v), B 2 ( V ) = 2 dfjH 2 (fj) 2(/(77)-l) 2 -3 + 3fi m (7?) C(r,,fj) , (49) 



" drj(e(fj) + l)c(r,,i}), -84(77) 



and 



df)C{r),ff) 



~,a 2 (f,), 



C( V ,fj)^g 2 (ij)a(fj) g( V )H(fj) - g (fj)^H( V ) 



a 2 (77)' 



with e(r?) = f 2 (rj)/n m (r]) and 



/(»?) = 1 + 



(50) 



(51) 



(52) 



g m denotes the value of g(r] m ), deep in the matter era before the cosmological constant was important. We also 
have a n i which denotes any primordial non-Gaussianity present. We shall set this to unity, representing a single field 
slow-roll inflationary model, in what follows. 
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We shall also require the perturbed energy density up to second-order 



^2 = 2^(2) _ 6i ^(2) _ 6jtf 2 $ (2) +24ff 2 $ 2 +6 ^2 + i6 $92$ 

a 2 

(where k 2 = 8irG), and the Laplacian of the perturbed velocity: 



+4tf(fi m - 2)$<9 2 $ - 4(Q m + 2)<j><9 2 $ - 4(3f2 m + 2)$9 2 $ - — $9 2 $ 



+ 



3a 2 H 2 



Hd 2 <& d z ^> + d 2 ^> d 2 ^> + Hd k <S> d 2 d k <S> + <9 fe $ d 2 d k ^> 



(53) 



(54) 



A. The averaged perturbed EFE 

By making use of the expansion at second order of the three dimensional volume element 



1 - 3* 



(V _ $ (2)^ 



(55) 



for any scalar function T, the Riemannian average (T)x> can be expanded in terms of the Euclidean average over the 
domain T> 



(T) 



h 
h 

Jv 



3 xT 

3 x 



(56) 



on the background space slices as: 

(T).p = T(°> + (T«) + (T< 2 >) + 3 [(T^X*) - (T* 1 )*) 



(57) 



where T' ), T^ 1 ' and T^ 2 ) denote respectively the background, first order and second order parts of the scalar function 
T = Tf(°) + TfW + Tf( 2 ). In the rest of the paper, every spatial average will be a Euclidean one. 
The averaged Hubble rate as defined by equation (25) is given by: 



H v = H- (*) 



2(l + z) s 

9H 2 n 

in. 



(ff(<9 2 $) + (<9 2 $)) + ($ 4) 



-3<$)<$) 



\2Hn m [jT<$ d 2 $) + ($ a 2 *)] + (1 + 3Cl m )H 2 (d k $ d k <S>) + (2 + 3fl m )H(d k ^ d k <i>) + (d fc $ a fe $)} 



2(1 + zf 

%H 2 n m 



i/($)(a 2 $) + ($)(9 2 $) 



-\{^) + \{l + z){d 2 vM). 
The averaged Friedmann equation (27) reads: 



(58) 



Hi = H 2 -\ (Q v -C v + n v ) - 2H{§) + 1(1 + z) 2 (d 2 <S>) - AH 2 (<S> 2 ) + 2H(<S> $) - \(l + z) 2 ($ d 2 <S>) 

Do o 



+ 



4(l + z) s 

9H 2 m 



(i + n m ) i? 2 (d fc $ + 2i7(a fe $ a fe $) + (a fc $ a fe $) -6#($)($) + 2(i + z) 2 ($)(<9 2 $) 



+H 2 (^) + K -(6 2 p), 
6 



(59) 



Finally, the averaged acceleration equation (28), which gives an effective Raychaudhuri equation, reads: 



3 aft*p 

a-D 



:m z 1 



+9H 2 (1 - Q m ) ($) + 3JT<$) - (1 + z) 2 (<9 2 $) 
+3i? 2 (9Q TO - 7) ($ 2 ) - 3tf ($ $) + (1 + z) 2 ($ <9 2 $) 
(1 + z) 2 



+ 



3H 2 fe 



+9ff + 27iJ 2 (1 - n m ) ($) 2 - 3(1 + z) 2 ($)(9 2 $) 

+3i/ 2 (i-|n m ) ($ (2) )-^(^>, 

The averaged curvature term is: 

K v = 2(l + z) 2 [2(9 2 $) + 6($a 2 $)+3(a fc $9 fe $) + 6($)(a 2 $) + (a 2 * (2) ) 

and the additional backreaction terms are: 
4(1 + z) 2 



H 2 (d k $ 9 fe $) + 2H(d k <S> d k <l>) + (d k <S> d k <i>) , 

3iJ($> + (1 + z) 2 (<9 2 $) - 15ff <$ $) - 3($ 2 > - (1 + zf [(* <9 2 $) + 2(<9 fc $ 8lfc$)] 
2(1 + z) 2 - 



3i? 2 Q r 



H(® d 2 <f>) + ($ <9 2 $) + 9ff ($)($) + 3(1 + z)^($)(<9 2 $) 



-i(i + z) 2 (a 2 a>( 2 )) + ^(<i>( 2 )), 



fCj) — 



8(1 + z) 2 

3i?^ m 

8(1 + z) 2 

3HQ m 
8(1 + z) 4 
27# 4 ^ 
8(1 + z) 2 

-2i/(i + z)(a 2 u (2) ) 

(1 + z) 2 (4iJ 



H(d 2 ^) + (a 2 $)j + 6($ 2 ) - 6($) 2 

2iJ($ <9 2 $) + 2($ <9 2 <j>) + 3H(d k <S> d k <S>) + 3(d k $ <9 fc $) 
H 2 (d 2 $) 2 + 2H(d 2 <f>)(d 2 <l>) + (a 2 $) 2 

-3i? 2 ($)(a 2 $) - 3ij($)(a 2 $) + ij($)(a 2 $) + ($)(a 2 $) 



H 2 tt 
2 

~3 

1 

+ 37T 
4(1 + z) 2 



h[i- jn m ) (d 2 <t>) + <s 2 <i>) 



h 2 (4 - 3n m ) ($ a 2 $) + 4ff($ a 2 $) + 3# ($ <9 2 $) + 3(6 a 2 $) 

3ff 2 (3^ - 2fl m - 4)(9 fe $ 9 fc $) - 8H(d k $> d k 4>) - 2(2 - 3n m )(d k <f> d k <f>) 
H 2 (d 2 $ <9 2 $) + 2H(d 2 <i> d 2 $>) + (9 2 $ <9 2 $) 



+ \h 2 (4 - m m ) ($) (<9 2 $) + 2(6) (o> 2 $) + 2iT($) (<9 2 $) + 4if ($) (d 2 $) 



+ 



4(l + z) s 

9ff 2 a„ 



H 2 (d 2 <S>} 2 + 2H(d 2 <S>)(d 2 <S>) + (a 2 *) 2 



+ + z){d 2 v {2) ) - ^H{1 + z)(c>V 2) ). 
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B. The averaging procedure 



We have expanded all perturbative quantities in terms of a Euclidean average, which will let us calculate the 
averaged quantities of interest in terms of the primordial power spectrum of Our Euclidean averaging procedure 
follows [6] which we summarise here. We apply our spatial average over a spherical domain using a finite window 
function, and then we apply an ensemble average to tell us what this spatial average will be in a typical domain. We 
may also calculate the variance we may expect as we move from domain to domain, which we would expect to vanish 
as the domain size becomes large. 

The spatial average involves a specific domain V which is arbitrary in principle. Here we use a window function W 
to specify the domain size and shape. For any spatial variable X(x), let 

(X(x)) = 1 J d 3 xW{x/R v )X{x) (66) 

where the integral is taken over all space, and Rx> specifies the size of the domain. We shall take our domain 
to be spherical as we expect a roughly isotropic distribution, and use a Gaussian window function for simplicity: 
W(y) = e-y 2 ' 2 , so that 

V = Jd 3 x W(x/R v ) = 4irRl J y 2 W(y)dy = (2 7 r) 3 / 2 J R|, (67) 

for any R. Here we use 

W{kR v ) = \ I d 3 xW(x/R v )e- lk x (68) 

as the Fourier transform of W(x/Rj>)/V. A top-hat window function could also be used instead. 

Without a specific realisation of perturbations to work from we must calculate what the average of X(x) would 
be in a typical domain of size Rv ■ Primordial perturbations are usually taken to be Gaussian fluctuations with zero 
mean, which implies that we can calculate what the "average of X(x) would be in a typical domain of size i?x>" by 
taking an ensemble average over many domains [38, 39]. Following [6] we denote this separate, additional average by 
an over-bar. 

It is clear that once the ensemble average is taken in our averaged backreaction equations, all stand-alone first-order 
perturbations will vanish if we assume that $(x) is a Gaussian variable, and only quadratic - second-order - quantities 
will remain. Let us therefore apply this spatial-then-ensemble average to the product of two Gaussian random scalars 
A(x) and B(x). We find first that 



(A(x)B(x)) = j^j d 3 ^i J d 3 fc 2 W(\k 1 +k 2 \Rv)A(k 1 )B(k 2 ). 



(69) 



Now, if A and B are statistically homogeneous, then we may write A{k) = A(k)E(k) where E is a unit random 
variable satisfying E*(k) = E(—k) and 

E(k 1 )E(k 2 ) = 5(k 1 + k 2 ). (70) 
Assuming that A and B are perfectly correlated random variables, we then find 

(A(x)B(x)) = ^jd 3 k A(k)B(k) . (71) 

Note that this is a constant as we'd expect, and that the window function has dropped out of the average. That is, 
once a statistical average is taken, the size of the domain doesn't make any difference for the product of Gaussian 
Random fields. We return to this below. Similarly we find 

(diAjx) WB(x)) = +^ j d 3 k k 2 A(k)B(k) . (72) 

We shall also require ensemble averages of the product of spatially averaged variables, which arise from commuting 
the spatial average with the time derivative. For example 

(A(x))(B(x)) = ^Jd 3 k W(kR v ) 2 A(k)B(k) , (73) 
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which does depend on the domain size - it is just (A(x) B{x)) with a squared window function stuck in the integral. 

Most of the terms we are dealing with are scalars schematically of the form d m Q(x)d n Q(x) where m and n 
represent the number of derivatives (not indices), such that m + n is even so that there arc no free indices. (For 
example, di$d 2 d l & has m = 1 and n = 3.) Then 



(d m <S>(x)d n <i>(x)) = (_i)(™+3«)/2 J dk fc m +"- 1 7' $ (fc). (74) 



In particular, we then also find 



(d 2 [d m $(x)d n $(x)]) = 0. (75) 



Time derivatives of the Bardeen potential may be dealt with using $(t, x) = — Q(t, x). 

1. Ensemble Average of Inverse Laplacian Terms 

In the solutions for the second-order Bardeen potentials, we encounter terms involving inverse Laplacians of 
quadratic first-order variables, such as 

d- 2 {di§d%) and 8r 2 d i di[Br i $d j §\. (76) 

These need some care when we take the ensemble average. 

We need to find the ensemble average of an object which is the inverse Laplacian of a product of Gaussian Random 
Fields. Consider 

8 2 X{x) = A(x)B(x) (77) 

so that we may write, formally, 

X(x) = d- 2 [A(x)B(x)}. (78) 

(Note that X(x) itself is not a product of Gaussian random fields since d 2 [A(x)B(x)] = - i.e., the Laplacian of the 
product of Gaussian random fields has zero mean - while d 2 X(x) ^ 0.) In Fourier space we have 

X{k) = (2^ fc ~ 2 / d " k ' A ^') B ^ k fe ')' ( 79 ) 

Taking the ensemble average of this equation gives us a 5(k)/k 2 term; if we are to find X(x) then we would have to 
integrate over fc-space which now involves this awkward divergent term. Let us start instead from the formal solution 
to Eq. (77) in real space which is, ignoring boundary terms and any homogeneous solution, 

- -SfSji Sw^r\ I ^ / ^e=.)B(* 2 )e«-*>»'. (81) 
Now, if A and B are Gaussian Random Fields, then we have 

1 1 f d 3 x' 



X ^ = -^(2^J W^\J d kA{k)m - m 
That is, the inverse Laplacian gives a divergent / d 3 x/x factor. To relate to Fourier space, we can use the relation 

1 1 r ik-(x-x') 

Sk e ——, (83) 



\x-x'\ 2tt 2 J k 

so that we have 



V - (2 ^ ):! / d 3 k f d*k l6 ^A{k')B[l')<' k \ (<S4) 



13 



which is what we get if we start in Fourier space. 3 (Despite appearances this isn't actually a function of x - we can 
expand out the k integral into angular and radial parts with the fc z -axis parallel to x to see this.) It is clear that 
calculating the ensemble average of an inverse Laplacian using Eq. (82) will be easier. 

Now we must find the ensemble averages of d^ 2 did : >[d t ^dj^], and d~ 4 did^[d i ^dj^ l \. Let us investigate 

cT" {did0A{x)&B{x)]} = J d 3 fcl J d 3 fc 2 A(fc 1 )i?(fc 2 - fci)/(fci,fc 2 )e ifca -" (85) 

where n — 2 or 4 and 

f(k 1 M) = {ki - k2){k ;- kl) ' k2 - (86) 

If we try to take the ensemble average at this stage we get 5(fc 2 )//c 2 as we ^ as mixed arguments in A and B (which 
means we can't easily perform the angular part of cither integral). We may however write S(|fe 2 — fc 1 |)<5(fc 2 ) = 
_B(fci)<5(fc 2 ) and perform the angular part of the fei-integral, which gives 4 



d-d i d j [^A(x)diB(x)] = J d 3 k'^e ik '* J™ dk k*A(k)B(k) 



(90) 



1 f°° 

- — ~ / dfc k 4 A(k)B(k) if n = 2 (91) 
/ 1 / dk k 4 A(k)B(k) if n = 4. (92) 

!7r) 3 7 |£C-£C'|7 



3(2 

Hence, for A = B = $, we find 

(d- 2 did 3 [d l ^(x)d^(x)}} = \ t dkkV^{k) = l(d i ^{x)d i ^{x)), (93) 

3 7 3 

Formally, for a scale-invariant spectrum the latter of these is infinite. However, where this appears in the second-order 
Bardeen potentials, it actually doesn't contribute, because we have shown the important result: 



9- 2 (frfbdi®) - 3d~ 4 did3 (d^dj®) = 0. (95) 



3 Using 5(k) = — *— / d A ye lk -y we find the identity / d 3 fc^ e lfe - = — f - d ^ / - 
(2-k) 6 J J k l 4ir J 



(2tt) 3 J J J k 2 An J \x + y\ 

4 This is a bit ambiguous, however, because /(fci, k 2 ) is singular at k 2 = (although doing this does give the correct answer, provided B 
is well behaved). For example, we could take some of / into this and write B(\k 2 — fei|)(fei • ^2)5(^2) = B(fci)(fei • 0)<5(fc2) = 0, which 
gives the wrong answer. To be more precise, let us write 

A(k) = J ' A a k'A(k~k')S(k') (87) 

so that we have 

d- n d i d j [d i A(x)d j B(x)] = ( ~ 2 ^" 3 /2 J d 3 fci J d 3 fc 2 j d 3 fc 3 A(fc! - k A )B(k 2 - k 1 )f(k 1 ,k 2 )S(k 3 )e lk ^ x (88) 
and now, taking the ensemble average, 

d-"-d i dj[d i A(x)djB( X )] = ^l"^ J d3fc i J d3fc 2 j d 3 fe 3 A(|fc 1 -fc3|)B(|fc2-fci|)/(fc 1 ,fc2)5(fc 3 )5(fc2-fc3)e lfe2 - a: 

= -J^^ J d3k 3j d 3 k 4 A(k 4 )B(k 4 )f(k 3 -k 4 ,k 3 )S(k 3 )e lh ^ (89) 

where the k 2 integral was done and k 4 = k 3 — k±. We now perform the angular part of the k 4 integral, giving the result stated in the 
text. 
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That is, if d 4 X(x) = d l d : > [d^(x)dj^(x)~\ then X(x) = (angle brackets on the indices denote symmetric trace 
free part of a tensor), which also says that the ensemble average of the scalar part of d^dj^ is zero. Hence, after 
ensemble averaging, the second-order Bardeen potentials have no divergent terms in the IR from inverse Laplacians. 
Similarly, using these results it is possible to show that 

(02tf(2)(a;)) = (0 2 *( 2 )(a!)) = (d2 v W(x)) = 0. (96) 



C. Length scales, smoothing and divergences 

There are two separate divergences which occur in the averaged equations. We have a UV divergence from terms 
like J dkk a V<j>, and IR divergences from the same types of terms. We discuss these in turn. 



1. UV divergences 

On small scales the transfer function behaves as (In k) / k 2 for pure CDM, and so integrals such as Eq. (74) diverge 
whenever m + n ^ 4. Recall that this is independent of the size of the domain which has already been taken into 
account. There arc different ways to deal with this divergence. If inflation introduces a tilt to the primordial power 
spectrum then the divergence only formally occurs for m+n > 4 - that is, in terms which don't appear here. However, 
in the case m + n = 4 for a tilt of n s = 0.95, for example, convergence of the integrals only happens for k/k eq > 10 20 
which is absurdly tiny scales - so this is effectively a divergence. Secondly, if the baryon fraction is high enough, this 
cuts off small-scale power and can remove some of the divergences. But we should expect finite results for basic CDM 
too (though of course in reality there is a cutoff at tiny galactic scales). Alternatively, we may introduce a small-scale 
cutoff by hand in the Fourier integrals. The corresponding sampling in real space is then oscillates a bit like a;sin(l/x) 
- which is clearly not suited for the purpose of averaging. 

Instead, we shall remove the UV divergences which appear by smoothing our perturbation &(x) in real space before 
averaging it. Since linear theory is only accurate up to some scale and under-estimates power below this scale, it makes 
sense to smooth away structure on scales smaller than this before we apply our averaging procedure - presumably 
higher-order perturbation theory will be required to correct for this. We replace <!> with a weighted average over 
nearby points using [39] 

*s(a:)= ^ Jd 3 x'W{\x'- x\Rs) *(x / ), (97) 

where for simplicity we use the same window function as we use in calculating the averages, but we separate the 
smoothing length scale R$ from the averaging length scale R V - In Fourier space this amounts to replacing 3>(fc) 
everywhere with W(kRs)$(k), and V$> with W(kRs) 2 V<s>. Note that the only terms which are made finite by this 
smoothing are the terms with four derivatives in them such as (d 2 Qd 2 Q) - although all terms are affected to some 
degree. These terms appear in the generalised Raychaudhuri equation, but not in the averaged Hubble rate. Hence, 
this divergence can significantly affect the effective equation of state, but not the averaged exp ansion rate. On the 
other hand, there are terms in both the Friedmann and Raychaudhuri equations like (9 2( f>) (d 2 Q) which have the same 
kind of divergence as Rx> — ► 0, but these divergences are controled by the domain size. 



2. IR divergences 

While the UV divergences are hopefully the result of using second-order perturbation theory and so can be removed 
by hand until a fully rclativistic renormaliscd approach is available, the IR divergences are maybe a much deeper 
issue (see, e.g., [40] for a discussion). For pure dust the integrals converge as k — > (the only divergent terms are 
(<I> 2 ) which cancel in the Friedmann equation); this IR divergence shows up in the variance of the expansion rate [6]. 
For ACDM, however, this doesn't happen because $ ^ and many terms appear which have an IR divergence. For 
a scale-invariant primordial power spectrum the IR divergence is only logarithmic, so is not really an issue. They 
may be removed by inserting a Hubble scale cutoff in the lower limits of the fc-integrals - but why choose the Hubble 
scale? 

During matter domination the Hubble or horizon scale is not really of direct physical significance, other than a 
symbolic 'size' of the universe. As far as structure formation goes it plays no role at all at this perturbative level. The 
only physical scale (crudely speaking) is the equality scale k eq which is the Hubble scale at matter-radiation equality. 
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The longest wavelength modes in the present-day universe are set by the first modes leaving the Hubble radius at 
the beginning of inflation, and not by the modes which happen to be entering the horizon today. Modes longer than 
the present day horizon are effectively homogeneous and so are usually considered to be part of the background, 
effectively rcnomalising them away. However, this argument only works if those super-Hubble wavelength modes, 
when renormaliscd to be part of the background, enter the field equations as curvature terms. This is not the case for 
backreaction. The time evolution of the averaged quantities is complicated and doesn't scale simply as I /a 2 . Hence, 
such modes may be considered part of the homogeneous background, but not in a theory of gravity obeying the usual 
Einstein field equations. Indeed, it is exactly that homogeneous background that the backreaction approach aims to 
calculate. 

It is not clear whether modes longer than the Hubble radius are real - if the standard inflationary picture is correct - 
and might have to be taken into account in the averaging backreaction approach. If the super-Hubble spectrum is 
tilted to the red this could be important [6] ; conceivably back-reaction may allow a mechanism to probe this part of 
the power spectrum. We set the lower IR limit of fc-integration to be ku/L where L is the largest wavelength mode 
we are prepared to include, in units of the present day Hubble scale. For simplicity we just assume a scale-invariant 
spectrum for the whole range of L, and fix L to be ten times larger than the Hubble scale. Generally we get the same 
answers as if we fix it at the Hubble scale, but this allows us to have a domain as large as the Hubble radius. 

IV. THE AVERAGED HUBBLE RATE AND DECELERATION PARAMETER 



We shall use length scales intrinsic to the model as reference points for smoothing and averaging: the Silk scale, 
-l 

cq > 



A-siik the equality scale, fc and the Hubble scale, fcjj 1 , and so the baryon fraction appears as this governs the Silk 
scale. Recall that [49] 



fc silk « i.6 (n h h 2 y- oz (n h 2 ) v "' [i + (io An h 2 y"- 95 \ m pc ~\ (98) 

k cq w 7.46 x KT^o^Mpc- 1 , and = JL M P C_1 > (99) 

ouuu 

where f2{, and flo are the baryon and total matter contributions today and Hq = 100 /ikms _1 Mpc _1 . 

We shall use two models for comparison: Einstein-de Sitter with h = 0.7 and 5% baryon fraction (WMAP5 [50] 
estimates w 0.046). This has k~ q ~ 27.9Mpc and ~ 6.0Mpc. The other model we shall use is a concordance 
model with fl = 0.26, h = 0.7, /baryon = 0.175 (this is the WMAP5 best fit [50]). The key length scales in this model 
are ~ 107.2Mpc and fc^ lk ~ 11.5Mpc. Both models have fc^ 1 ~ 4.3Gpc. To calculate the integrals we use transfer 
functions presented in [49]. All lengths shown are in Mpc. 

We set L — 10; that is, all fc-integrals have an IR cut-off set at ten times the Hubble scale. 



A. Hubble rate 



There are different aspects of the backreaction we wish to probe, and some subtleties arise because we have to 
take an ensemble average of our equations. When we examine the Hubble rate we are interested in two things: the 
dynamics of the expansion rate, and the averaged Friedmann equation. In the Friedmann equation we are interested 
in quantifying the new terms which enter the field equations as a result of averaging, which are the new components 
which drive the expansion; within the context of dark energy, it is common to think of these as effective fluid or 
curvature terms. While the spatial average of these two agree up to perturbative order, when we take the ensemble 
average we ascertain different information. 

Take the ensemble average of the Hubble rate given by the generalised Friedmann equation. For a given domain 
size this tells us the expectation value of the averaged Hubble rate we might expect to find dynamically in the field 
equations. When we present Hz below, we have usually calculated 




i.e., we have taken the ensemble average of the Friedmann equation and then taken the square-root. This does not 
yield the same answer as taking Hz using Eq. (58) directly, which is the expectation value of the kinematical Hubble 
rate (but note that if we square Eq (58), take the ensemble average, and then take the square-root we get the Hubble 
rate as calculated directly from the Friedmann equation). The difference of course is the 'ensemble- variance' of the 



1G 

Hubble rate, which may be defined by 

Y&x[H v ]= Hi - H% . (101) 

[Eq. (59)] [Eq. (58)] 2 



When we write [Eq. (58)] , this is developed to the correct perturbative order. 




i 1 1 1 1 1 1 1 1 1 r 



1 2 3 4 5 

z 

FIG. 1: The averaged Hubble rate as a function of redshift as a fractional change to the background Hubble rate. The top two 
curves show the Hubble rate Ht> calculated from the ensemble-averaged Friedmann equation, and the bottom two show the 
Hubble rate calculated directly, H x>- Both concordance and EdS models are considered, and the domain size is set at the Silk 
scale, and the smoothing scale is set to zero. 
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FIG. 2: Plots for Ht> with Rt> = l/fc eq uaiity, Rs = 0, with the variance included, as a function of redshift. 
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FIG. 3: The averaged Hubble rate today as a function of domain size, with the smoothing scale removed, Rs = 0, and with 
the variance included. 
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FIG. 4: The averaged Hubble rates today as a function of domain size, and the variance (top), with the smoothing scale removed, 
Rs = 0. Note that both EdS and concordance models become independent of domain size when it becomes sufficiently large. 
Intriguingly this happens around 100 — 200Mpc; at this scale, for both models, we also have Hx> ~ Hx> and the variance reduces 
to the fiducial 10 -5 level. 



18 



250 



200 



150- 



100- 



50 




scale at which Hp, becomes roughly 
independent of R-p and H-p ~ Hp> 



scale when '2cr' uncertainty in H-p is < 1% of Hq 



0.1 0^2 03 o'.4 ' o'.5 o'.6 07 0^8 o'.9 l'.O 



equality scale, 1/A; equa i it y 
Silk scale, l/fcsilk 



FIG. 5: Homogeneity scale in the averaged Hubble rate, as a function of present-day matter density. As the domain size is 
increased the averaged Hubble rate settles down to a constant value which is different from Ha (see Fig. 4). We may investigate 
this scale as a function of fio by deciding when Hx> is sufficiently independent io Rt> (here, \d\n(H-D / Hq — l)/dlni?x)| < e - 
i.e., when the slope in the log-log plot of Fig. 4 levels off). To illustrate, we have arbitrarily fixed this e so that for fio = 0.1, 
the homogeneity scale coincides with the equality scale, but this curve can be moved up and down to taste. Finally, another 
useful scale is when the variance becomes negligible: here we show the scale for which 2-i/varpfn] = 0.01i/o- Note how the 
latter of these measures behaves in the opposite way to the scales intrinsic to the models - the equality and Silk scales. We 
have fixed the baryon fraction at 5% here. 



In Figs. 1 - 5 we show various aspects of the averaged Hubble rates. The backreaction effect grows during dust 
domination and decays after the dark energy transition. In Einstcin-de Sitter, then, the backreaction effect is largest 
today, while in the concordance cosmology it peaks around z ~ 0.5. The smoothing scale and most importantly the 
averaging domain size are crucial in deciding how large the effect of backreaction is. We sec also that whether we 
consider Hx> or Hp, matters considerably if the domain size is small; only if it is larger than a few hundred Mpc do 
these agree. 

As far as the Hubble rate is concerned there is no UV divergence in any of the integrals and so we have, for 
simplicity, set the smoothing scale to zero. In Figures 3 and 4 we show the dependence on domain size on the Hubble 
rate today. It is largest for small domain size, but we see that it levels off at a constant value for R-p bigger than 
a hundred Mpc, or thereabouts, and both Hubble rates are the same (i.e., the variance is sufficiently small). Recall 
that the domain scale only affects terms of the form (• ••)(•••), so as R-p becomes large those terms disappear, leaving 
terms like (• • •) which are independent of the domain scale. 

This can be interpreted by saying that there exists a scale which corresponds to the scale of homogeneity of the 
distribution of matter: there's an upper limit to domain size beyond which variables are fixed. This is explored 
further in Fig. 5, where we also show the scale beyond which uncertainty in H drops to negligible levels. The scale 
of homogeneity that comes out of the averaging process is similar to the equality scale, but behaves differently with 
fioj so the relation is not direct. In magnitude, this is similar to the scale of homogeneity inferred from large scale 
structure (LSS) surveys like SDSS [44, 45], but there are conflicting results in this area; see [46-48], where it is found 
that homogeneity cannot be inferred from LSS surveys out to 100/i _1 Mpc. We are not able to say in this study 
whether back-reaction effects could account for such a discrepancy, but it is possible that the naive scale derived in 
linear theory is refined significantly when non-linear backreaction effects are properly accounted for. This deserves 
further investigation. 

Thus we see that above a certain scale the backreaction leaves a small, scale-invariant, residue. Above the scale of 
homogeneity, the backreaction effect becomes scale independent, but not at a zero value; in other words, the effective 
homogeneous model of the perturbed Universe that is obtained by averaging on scales bigger than a few hundred Mpc 
exhibits a 'renormalized' background with non-zero backreaction. This is an important fact since we started with a 
given background, but we don't recover this background after averaging; that essentially means that stricto sensu, 
backreaction cannot be said to vanish on very large scales. In principle, this could provide a novel way of constraining 
the matter density by using measurements of this scale of homogeneity in large scale structure data. 
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Our results are quantitatively comparable to those of [18] obtained in the synchronous gauge. In particular, the 
variance effect due to the backreaction exhibits the same scale behavior, with comparable amplitude. Our estimate 
of the variance is consistent with the variance in the Hubble rate calculated directly from the variance of peculiar 
velocities in [41]. This is important for future research aiming at 1% uncertainty in measuring H [42, 43]. We find 
comparable results for the change to the Hubble rate in Refs. [20, 21] where it was found that the backreaction in 
the Poisson gauge gives a change to the Hubble rate today <~ 1CP 5 ; we find about the same for the concordance 
model when averaging on Hubble scales, but we have shown that the backreaction can be large when averaging on 
sub-Hubble scales. Any difference can be attributed in part to the different definition used to calculate the efffective 
Hubble rate, and we also employ a different averaging scheme. We also include the second-order Bardeen potentials 
in our analysis, which were neglected in [20, 21] though these are a sub-dominant contribution to the backreaction. 



B. The Raychaudhuri equation and deceleration parameter 

While we have seen the change due to backreaction in the Hubble rate arc small, and, in particular, have no UV 
divergence, we can get further information about the effect of the backreaction from perturbations by looking at the 
Raychaudhuri equation. These are sourced in part by the new backreaction terms such as K-x> and so on, but with 
many other quantities contributing. In terms of domain and smoothing scale, JCj> diverges as Rs — ► 0, and Qx> — C-d 
diverges as i?x> — > 0; it becomes domain size-invariant beyond a few hundred Mpc. All the other backreaction terms 
are fairly independent of domain and smoothing scales. 

We define an averaged deceleration parameter as 

M*) = -7^ (102) 

where a^/a-D is given by the generalised Raychaudhuri equation, above. To calculate the ensemble average, we first 
note that 



^ = H{zf ( 1 - ln m (z)) + SWR + 5™R (103) 
a-v \ 2 J 

where S^R are all first-order terms such as (<&) (which have zero ensemble average), S^R contains all second-order 
terms such as (<I> 2 ) and so on; similarly for the Friedmann equation = H 2 + S^F + S^F, where S^F = 0. We 
then find 



3^ , N lap,/, 3^ , N \ Hi SWRSWF 



= -1 + ^ m (z) - — ^ + ^1 - -Q m (z) j ^ + — 4 . (104) 

(To high accuracy we get the same result if we just replace with H v in Eq. (102).) To calculate the variance of 
qx> we must expand g|, to second-order to find q 2 -,. We then find 

Var[<ft>] - H- 4 1 1 [2 - 3n m (z)} 2 SW F 2 - [2 - 3Q m (z)]SW FS^ R + SWR 2 ^ . (105) 
In Figs. 6 and 7 we explore the behaviour of the expectation value of the deceleration parameter, and its variance. 



V. CONCLUSION 



In this paper, we have calculated the effect of the backreaction arising from averaging a perturbed FLRW spacetime 
in the Poisson gauge consistenly up to second order in perturbation theory. Inclusion of a cosmological constant gives 
the analysis additional intricacy as we include the second-order Bardeen potentials, complications which have been 
ignored in previous studies. In particular they contain non-local terms which require some care to show that they 
vanish once an ensemble average is taken. We utilised an averaging procedure which has two elements: a spatial 
average over a typical domain, and an ensemble average over domains. The ensemble average gives the expectation 
value for a quantity averaged on a domain V, and the variance gives an intrinsic uncertainty to the expected value. 
We discussed in detail the changes that arise from averaging the Hubble parameter and the Friedmann equation, and 
the deceleration parameter. 

Our key results are: 
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FIG. 6: The ensemble-averaged deceleration parameter and its variance, for the EdS and concordance models. We can see from 
the graphs on the left that the difference from normal behaviour grows during dust domination, and decays when the dark 
energy accelerated the expansion rate. Interestingly, the variance also starts to decrease then, reaching a minimum for z fa 0.22. 
Note the divergence with the smoothing scale becomes irrelevant if it is bigger than tens of Mpc (right). The dependence on 
domain scale depends crucially on the model, however (middle), and no longer matters if Rd < Rs (shown for the concordance 
model only, where the curves flatten off). 



• There exists a homogeneity scale inherent in the averaging procedure which is not simply the equality scale 
for a general f^o- As the domain size increases the averaged Hubble rate becomes independent of averaging 
scale, and the variance of the Hubble rate becomes negligible. We also see a similar scale in the behaviour 
of the deceleration parameter. While this scale is comparable to the one inferred from large-scale structure 
surveys [44, 45], we have seen that backreaction can give the homogeneity scale different to the equality scale 
somewhat. Further work could look into whether the homogeneity scale is calculated sufficiently accurately in 
linear theory, as backreaction may account for discrepancies in determining it from observations [46-48] . 

• For domains larger than the homogeneity scale, the backreaction does not go to zero; rather, there is a small 
residue left from the averaging, of order 10 -5 for averaging the Hubble expansion on Hubble-scale patches. 
Thus, the background is renormalised by structure formation, and so the 'background' at late times is not the 
same as the background at the end of inflation. 

• It is not clear exactly how large the change to the deceleration parameter is. A seemingly sensible and conser- 
vative choice of smoothing scale - the Silk scale - gives a sizeable change to the deceleration parameter, of the 
order of 10% (for our concordance model - more as f2o increases) , even when averaging above the homogeneity 
scale. If it is this large, this effect could be critical in determining the nature of dark energy. Furthermore, the 
deceleration parameter is pushed to more positive values by backreaction. The existence of a UV divergence 
in the effective Raychaudhuri equation implies that to quantify this properly, higher-order perturbation theory 
might be required. However, it has been argued [51] that going to higher-order in perturbation theory may 
make the UV divergence worse, not better. 5 



5 The problem is that at higher orders in perturbation theory higher derivatives of the first-order Bardeen potential appear, which, on 
the face of it, may make any UV divergences worse; it may even be argued that perturbation theory is not convergent on the basis of 
this [51, 52]. Our models of the universe may not be analytic. Alternatively, it may be more like convergence in the power series of 
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Rs 

FIG. 7: The ensemble-averaged deceleration parameter today as a function of smoothing scale (in Mpc). This gives the residual 
contribution to the deceleration parameter for very large domains - clearly the smoothing scale is crucial, and suggests that 
higher-oder contributions might be necessary in order quantify the backreaction effect in go- 



• Backreaction also introduces an uncertainty in the value of the Hubble rate and deceleration parameter which 
manifests itself through a variance in the value in different averaged regions. This variance has been estimated 
here and is found to be quite large if the averaging scale is below the homogeneity scale discussed above (around 
a few percent for small averaging scales), in approximate agreement with estimations previously made in the 
synchronous gauge [18], as well as estimates made independently of the relativistic averaging scheme [41]. This 
may be an important consideration for observationally determining cosmological parameters such as Hq to 
accuracy higher than a few percent. 

There exists an IR divergence inherent in this analysis. This means that all the modes bigger than the one 
corresponding to the averaging scale contribute to the backreaction effect, irrespective to their size, for a scale- 
invariant spectrum. In particular, the Hubble radius does not appear as a natural cut-off, and the smallness of 
the effect noted above only comes from the introduction of an arbitrary cut-off: by pushing this cut-off towards 
higher values, the backreaction effect is larger, though only logarithmically with scale. For a red spectrum, however, 
this divergence is power-law so could be important. In contrast to [6] we find that this divergence appears in the 
backreaction effect directly and not just the variance. The longest wavelength mode should perhaps be set by the 
start of inflation, though it is by no means clear exactly how to implement this properly. 

The assumption of a primordial spectrum of adiabatic Gaussian perturbations was also critical in this analysis, 
and all conclusions we present are reliant on this fact. In particular, because the first-order contribution to the 
backreaction vanishes for Gaussian perturbations, non-Gaussian initial perturbations may produce a much larger 
backreaction effect, of the order of the square-root of the variance we have presented here. This is an important 
outstanding issue to explore. 
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